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Abstract 

A version of the Dynamical Systems Method (DSM) for solving ill-conditioned linear 
algebraic systems is studied in this paper. An a priori and a posteriori stopping rules 
are justified. An algorithm for computing the solution using a spectral decomposition 
of the left-hand side matrix is proposed. Numerical results show that when a spectral 
decompositon of the left-hand side matrix is available or not computationally expensive 
to obtain the new method can be considered as an alternative to the Variational 
Regularization. 
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1 Introduction 

The Dynamical Systems Method (DSM) was systematically introduced and investigated 
in [11] as a general method for solving operator equations, linear and nonlinear, especially 
ill-posed operator equations. In several recent publications various versions of the DSM, 
proposed in were shown to be as efficient and economical as variational regularization 
methods. This was demonstrated, for example, for the problems of solving ill-conditioned 
linear algebraic systems [2] , and stable numerical differentiation of noisy data [8] , [9] , [3] . 

The aim of this paper is to formulate a version of the DSM gradient method for solving 
ill-posed linear equations and to demonstrate numerical efficiency of this method. There 
is a large literature on iterative regularization methods. These methods can be derived 
from a suitable version of the DSM by a discretization (see [H]). In the Gauss-Newton- 
type version of the DSM one has to invert some linear operator, which is an expensive 
procedure. The same is true for regularized Newton-type versions of the DSM and of their 
iterative counterparts. In contrast, the DSM gradient method we study in this paper does 
not require inversion of operators. 
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We want to solve equation 

Au = /, (1) 

where A is a linear bounded operator in a Hilbert space H. We assume that ([1]) has a 
solution, possibly nonunique, and denote by y the unique minimal-norm solution to ([T]), 
y 1 Af := M{A) := {u : Au = 0}, Ay = f. We assume that the range of A, R(A), 
is not closed, so problem ([I]) is ill-posed. Let f$, \\f — fg\\ < 5, be the noisy data. We 
want to construct a stable approximation of y, given {5, f$,A}. There are many methods 
for doing this, see, e.g., [4], [5], [6], [TT], [13], to mention a few books, where variational 
regularization, quasisolutions, quasiinversion, iterative regularization, and the DSM are 
studied. 

The DSM version we study in this paper consists of solving the Cauchy problem 

ObU 

u(t) = -A*(Au(t) - /), «(0)=uo, u ±N, u:=—, (2) 

where A* is the adjoint to operator A, and proving the existence of the limit lim^oo u(t) = 
u(oo), and the relation u(oo) = y, i.e., 

lim \\u(t) - y\\ = 0. (3) 

t— >oo 

If the noisy data f$ are given, then we solve the problem 

u s (t) = -A*(Au s (t) - f s ), u s (0) = u , (4) 
and prove that, for a suitable stopping time ts, and us := u$(ts), one has 

lim \\u$ — y\\ = 0. (5) 

In Section 2 these results are formulated precisely and recipes for choosing ts are 
proposed. 

The novel results in this paper include the proof of the discrepancy principle (Theorem 
3), an efficient method for computing us(t$) (Section [3|), and an a priori stopping rule 
(Theorem 2). 

Our presentation is essentially self-contained. 



2 Results 

Suppose A : H — > H is a linear bounded operator in a Hilbert space H. Assume that 
equation 

Au = f (6) 

has a solution not necessarily unique. Denote by y the unique minimal- norm solution i.e., 
y _L N := Af(A). Consider the following Dynamical Systems Method (DSM) 

u = -A*(Au-f), 

u(0) = U Q , 
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where uq _L M is arbitrary. Denote T := A* A, Q := AA*. The unique solution to © is 

u(t) = e- tT u + e- tT f e sT dsA*f. 

Jo 

Let us show that any ill-posed linear equation © with exact data can be solved by the 
DSM. 

2.1 Exact data 

Theorem 1 Suppose uq _L N . Then problem (|7|) has a unique solution defined on [0, oo), 
and u(oo) = y, where u(oo) = lim^oo u(t). 

Proof. Denote w := u(t) — y, wq = w(0). Note that wq _L J\f. One has 

w = -Tw, T = A* A. (8) 
The unique solution to (JSj) is w = e~ tT uiQ. Thus, 

\\w\\ 2 = /" " e- 2tx d{E x w ,w ). 
Jo 

where (u, v) is the inner product in H, and E\ is the resolution of the identity of the 
selfadjoint operator T. Thus, 

r\\n 

\\w(oo)f = lim / e~ 2tx d(E x w ,w ) = ||P^ || 2 = 0, 
where Pj^ = Eq — E_o is the orthogonal projector onto M. Theorem Q] is proved. □ 

2.2 Noisy data f s 

Let us solve stably equation ([6]) assuming that / is not known, but fg, the noisy data, are 
known, where ||/<5 — /|| < 6. Consider the following DSM 

u s = -A*(Au s - fs), u s (0) = u . 

Denote 

w s :=u s -y, T := A* A, w s {0) = w := u - y G M 1 . 
Let us prove the following result: 

Theorem 2 If lim^o ts = oo, lim^o^ = 0, and wq _L N , then 

lim 11^(^)11 = 0. 
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Proof. One has 

w s = -Tw s +r)s, rjs = A*(f 5 - /), H^H < \\A\\S. 
The unique solution of equation ([9|) is 

w s (t) = e- tT w s (0) + f e~^ T Vs ds. 
Jo 

Let us show that lim^oo = 0. One has 



lim ||u>5(t)|| < lim \\e tT w 5 (0)\\ + lim 

t— >oo t— >oo 4— »oo 



One uses the spectral theorem and gets: 



<t-s)T d 



f I" dE xm e-^ x ds 
Jo Jo 

dExm = 



svs 



\T\\ p t\ _ i 



mi ! _ e -u 



dE x r] S - 



Note that 



1 - e"' A 

< < t, VA > 0, t > 0, 

A 



since 1 — x < e x for x > 0. From (jlip and (fT2]h one obtains 



l-e-* A |2 



d(E x Vs,ri5) 



imi 



d{Exr]s,Vs) 



,2 II i|2 



Since 1 1 775 1 1 < ||-A||<5, from (fTUl) and (fl3j) . one gets 



\im\\w s (t 5 )\\ < lim ( \\e- t6T w 5 (0)\\ +t s 6\\A\ 



Here we have used the relation: 



lim \\e- tsl w s (0)\\ = \\Pjs/w \\=0, 

and the last equality holds because wq G Af ± . Theorem [2] is proved. 

From Theorem [21 it follows that the relation tg = 7 = const, 7 € (0, 1) and 
is a constant, can be used as an a priori stopping rule, i.e., for such t$ one has 

lim \\u s (t s ) - y\\ = 0. 

o— >0 



2.3 Discrepancy principle 

Let us consider equation © with noisy data fs, and a DSM of the form 

us = -A*Au s + A*f s , u s {0) = u . (15) 

for solving this equation. Equation (fT5|) has been used in Section [221 Recall that y denotes 
the minimal-norm solution of equation ©. 

Theorem 3 Assume that \\Auq — f$\\ > CS. The solution tg to the equation 

h(t) := \\Au s (t) -fs\\ = CS, C = const, C G (1, 2), (16) 
does exist, is unique, and 

lim \\u 5 {t s )-y\\ =0. (17) 

Proof. Denote 

v s (t) := Au s (t) - f 5 , T:=A*A, Q = AA*, w{t) := u(t) - y, w :=u -y. 
One has 

iMt)|| 2 = 2Re(Au s (t), Au 5 (t) - f s ) 
dt 

= 2Re(A[-A*(Au s (t) - fs)],Au s (t) - f 5 ) ^ 
= -2\\A*v s (t)\\ 2 <0. 



Thus, ||V(s(i)|| is a nonincreasing function. Let us prove that equation (|16f) has a solution 
for C G (1,2). Recall the known commutation formulas: 

e~ sT A* = A*e~ sQ , Ae~ sT = e~ tQ A. 

Using these formulas and the representation 

u s (t) = e- tT u + f e-^ T A*f s ds, 
Jo 



one gets: 



v s (t) = Au s (t) - fs 

= Ae- tT u + A [ e~^ T A*fsds - f 5 
Jo 

= e- tQ Au + e~ tQ [ e sQ dsQf s - fs 
Jo 

= e- tQ A(u -y) + e~ tQ f + e' tQ (e tQ - I)f s - fs 
= e' tQ Aw - e~ tQ fs + e' tQ f. 
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Note that 

lim e~ tQ Aw = lim Ae~ tT w = APj^wq = 0. 



Here the continuity of A, and the following relation 

IIT|| 

^ im e~ tT wo = lim / 



r\\T\\ 

lim e~ tT wo = lim / e~ st dE s wo = (Eq — E-q)wo = Prfwo, 



were used. Therefore, 

lim \\v s (t)\\ = lim ||e- t( 3(/ - Mil < ||/ - Ml < 5, (19) 

i— >oo t— >oo 

because ||e _t< 5|| < 1. The function h(t) is continuous on [0,oo), h(0) = \\Auq — Ml > CS, 
h(oo) < 5. Thus, equation ([TBI) must have a solution t$. 

Let us prove the uniqueness of t$. Without loss of generality we can assume that 
there exists t\ > tg such that HAu^ti) — fs\\ = CS. Since ||i>5(£)|| is nonincreasing and 
||u 5 (t 5 )|| = ||v 5 (ti)||, one has 

||«i(t)|| = \\v 5 (ts)\\, Vie [t s ,h]. 

Thus, 

^IK(£)|| 2 = o, vte(t{,*i). (20) 

at 

Using (fTSj) and ([2D]) one obtains 

A*vg(t) = A*(Au s (t) - f s ) = 0, Vi G [ti,ti]. 

This and (fT5|) imply 



u 5 (i)=0, ViG(t 5 ,ii). (21) 



One has 



ug(t) = -Tu s (t) + A*f 5 

= -T^e- tT u + J e~^ T A*f 5 ds^j + A*f s 

= -Te- tT u -(I- e' tT )A*fg + A* f s 
= -e- tT (Tu -A*f s ). 



(22) 



From flSSJ) and flST}, one gets Tu - A*/ = e ff e _ff (Tu - A*/) = 0. Note that the 
operator e tT is an isomorphism for any fixed t since T is selfadjoint and bounded. Since 
Tu - A*f = 0, by (JS2J) one has u 5 (i) = 0, u 5 (i) = it a (0), Vi > 0. Consequently, 

C<5 < p««(0) - M| = \\Aus(t s ) - Ml = CS. 

This is a contradiction which proves the uniqueness of tg. 
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Let us prove (|17|) . First, we have the following estimate: 

\\Au(t s ) - /|| < \\Au(t s ) - Au 5 (ts)\\ + \\Au s (t s ) - f s \\ + \\f s - f\\ 



< 



Q I" e sQ Qds 
Jo 



\\fs-f\\ + C5 + 6. 



(23) 



Let us use the inequality: 



-t s Q I S ^QQds\\ = \\I - e ~ tsQ \\ < 2, 
JO 



and conclude from (I23H. that 



Secondly, we claim that 



hm\\Au(t 5 )-f\\ = 0. (24) 

lim tx = oo. 

<5^0 

Assume the contrary. Then there exist to > and a sequence (i<5 n )J£Li, £<s n < to, such that 

lim \\Au(t 5n )-f\\=0. (25) 

n— >oo 

Analogously to (fTBj) . one proves that 

cZ|M' 2 



eft 



< 0, 



where := Au(t) — f. Thus, ||u(t)|| is nonincreasing. This and (|25p imply the relation 
|Kt )|| = ||M*o)-/||=0. Thus, 

= v(t ) = e- toQ A(u -y). 

This implies A(u - y) = e toQ e~ toQ A(u - y) = 0, so u - y £ M. Since u - y £ Af- 1 -, it 
follows that uq = y. This is a contradiction because 

C5 < \\Au - f s \\ = \\f - f s \\ < 6, 1<C<2. 

Thus, lim^o^a = oo. 

Let us continue the proof of (fT7|) . Let := us(t) — y. We claim that ||iu,5(i)|| is 

nonincreasing on [0, tg]. One has 

^IN(*)I| 2 = 2Re{u s (t),u s (t)-y) 

= 2Re(- A* (Au 5 (t) - f 5 ),u 5 (t) - y) 

= -2Re(Au 5 {t) - fs, Au 5 (t) - fs + fs - Ay) 

< -2\\Au 5 (t) - f s \\ (\\Au s (t) -f 5 \\- \\f s - f\\ 



< 0. 
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Here we have used the inequalities: 



\\Au s (t) - f s \\ >C5> \\f s - Ay\\ = 5, Vi 6 [0,t s ]. 

Let e > be arbitrary small. Since lim^oo u(t) = y, there exists to > 0, independent 
of 6, such that 

IK*o)-y|| < |- (26) 

Since lining ts = °°, there exists So such that tg > to, \/5 G (0,5q). Since ||u>,5(t)|| is 
nonincr easing on [0, tg] one has 

\\w s (tg)\\ < \\wg(t )\\ < ||ii«(to) - «(to)|| + ||t*(to) - Z/H, V5G (0,<Jo). (27) 

Note that 

IM*o) - «(to)|| = Ik"* 07 f° e sT dsA*(f 5 - f)\\ < \\ e - toT f° e sT dsA*\\5. (28) 

JO JO 

bmce e 

~ toT I t0 e sT dsA* is a bounded operator for any fixed to, one concludes from ([25]) 
that lim^o ll^^o) — ' u (*o)|| = 0- Hence, there exists 5\ £ (0, So) such that 

\\ug(t ) -u(t )\\ < |, V<5 6(0,<5i). (29) 

From (I2BD-(I2"91. one obtains 



||««(*«) " V\\ = \\Mts)\\ < \ + \ = e, G (0,<5i). 
This means that lim^o u s(ts) = V- Theorem [3] is proved. □ 



3 Computing us(ts) 

3.1 Systems with known spectral decomposition 

One way to solve the Cauchy problem (|15p is to use explicit Euler or Runge-Kutta methods 
with a constant or adaptive stepsize h. However, stepsize h for solving (|15p by explicit 
numerical methods is often smaller than 1 and the stopping time tg = nh may be large. 
Therefore, the computation time, characterized by the number of iterations n, for this 
approach may be large. This fact is also reported in [2J, where one of the most efficient 
numerical methods for solving ordinary differential equations (ODEs), the DOPRI45 (see 
PQ), is used for solving a Cauchy problem in a DSM. Indeed, the use of explicit Euler 
method leads to a Landweber iteration which is known for slow convergence. Thus, it 
may be computationally expensive to compute ug(tg) by numerical methods for ODEs. 

However, when A in (|15j) is a matrix and a decomposition A = USV* , where U and V 
are unitary matrices and S is a diagonal matrix, is known, it is possible to compute ug(tg) 
at a speed comparable to other methods such as the variational regularization (VR) as it 
will be shown below. 
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We have ^ 

u s {t) = e- tT u + e- tT [ e sT dsA*f s , T := A* A. (30) 
Jo 

Suppose that a decomposition 

A = USV*, (31) 

where U and V are unitary matrices and S is a diagonal matrix is known. These matrices 
possibly contain complex entries. Thus, T = A* A = VSSV* and e T = e vssv * . Using the 
formula e vssv * = Ve ss V*, which is valid if V is unitary and SS is diagonal, equation (|30p 
can be rewritten as 

u s (t) = Ve- tSs V*u + V f e {s - t)Ss dsSU*f s . (32) 

Jo 

Here, the overbar stands for complex conjugation. Choose uq = 0. Then 

u s (t) = V f e^ Ss dsSh s , h 5 := U*f s . (33) 
Jo 

Let us assume that 

A*f s + 0. (34) 



This is a natural assumption. Indeed, if A* f$ = 0, then by the definition of hg in (|33p . 
relation V*V = I, and equation ([3T]) . one gets 

Sh s = SU*f 5 = V*VSU*f s = V*A*f s = 0. (35) 

Equations ([35]) and (f33j) imply Ug(t) = 0. 

The stopping time tg we choose by the following discrepancy principle: 



\\Au 5 (ts) - f 5 \\ 



5 e {s - ts)Ss dsSSh 5 - h s 



-^^^H = CS. 



where 1 < C < 2. 

Let us find tg from the equation 



<f>(t) := i/>(t) - C6 = 0, ip(t) :=\\e- tbb hg\\. (36) 

The existence and uniqueness of the solution tg to equation (|36p follow from Theorem [3l 
We claim that equation (1361) can 6e solved by using Newton's iteration ([43]) for any 
initial value to such that 4>(to) > 0. 

Let us prove this claim. It is sufficient to prove that <j)(t) is a monotone strictly convex 
function. This is proved below. 

Without loss of generality, we can assume that hs (see (j36j) ) is a vector with real 
components. The proof remained essentially the same for hg with complex components. 
First, we claim that 

tSSi 



SShg^O, and \\V SSe~ tbb h 5 \\ + 0, (37) 
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so ij)(t) > 0. 

Indeed, since e~ tss is an isomorphism and e~ tss commutes with V SS one concludes 
that \\y/SSe- t§s hs\\ = iff VSSh 5 = 0. If VSSh s = then Sh s = 0, and then, by 
equation ([35]) . ^4*/5 = Shs = 0. This contradicts to the assumption (j34l) . 

Let us now prove that (j) monotonically decays and is strictly convex. Then our claim 
will be proved. 

One has 

±(e- t§s h s ,e- t§s h s ) = -2(e- t§s h s ,SSe- tBs h s ). 
at 

Thus, 

Equation ([MI), relation ((37]), and the fact that (e- tSs h s , SSe~ tSs hs) = \\VsS e - tSs h 5 \\ 2 
imply 

ip(t) < 0. (39) 



From equation ([38]) and the definition of ip in ([36]) . one gets 

= -(e-* 55 ^, S^e"' 55 /^) (40) 

Differentiating equation (|40|) with respect to t, one obtains 

V?(t)#) +V> 2 (i) = (SSe' tSs h 5 ,SSe- tSs h 5 ) + (e- tSs h s ,SSSSe- t§s h s ) 
= 2\\SSe- tSs h s \\ 2 . 

This equation and equation (|38p imply 

l p -tSS h: QQ p -tSS h \2 

mm = 2\\sse- iss h 8 f - ( ,: t f Sh „ 2 f) > n^-^^n 2 > o. (4i) 

|| 6 <i(5|| 

Here the inequality: (e~ tss hs,SSe~ tss h s ) < \\e~ tss hs\\\\SSe~ tss hs\\ was used. Since 
tp > 0, inequality ([IT]) implies 

> 0. (42) 

It follows from inequalities ([39]) and ([4"2]) that 0(i) is a strictly convex and decreasing 
function on (0, oo). Therefore, tg can be found by Newton's iterations: 

. _ . 4>(tn) 



*(*») 3 

~ n+ (SSe-^S hs;e -t n SS hs) W e ^ n-U,l,..., 

for any initial guess to of tg such that (j){to) > 0. Once ^ is found, the solution ug(tg) is 
computed by ([33]) . 
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Remark 1 In the decomposition A = VSU* we do not assume that U, V and S are 
matrices with real entries. The singular value decomposition (SVD) is a particular case of 
this decomposition. 

It is computationally expensive to get the SVD of a matrix in general. However, there 
are many problems in which the decomposition (131D can be computed fast using the fast 
Fourier transform (FFT). Examples include image restoration problems with circulant 
block matrices (see [7]) and deconvolution problems, (see Section [42]) . 

3.2 On the choice of t 

Let us discuss a strategy for choosing the initial value to in Newton's iterations for finding 
tg. We choose to satisfying condition: 



by the following strategy 

1. Choose to := as an initial guess for to- 

2. Compute <f)(to). If to satisfying ([34"|) we are done. Otherwise, we go to step 3. 

3. If 4>{to) < and the inequality (f>(to) > 5 has not occurred in iteration, we replace to 
by |^ and go back to step 2. If 4>{to) < and the inequality 4>{to) > 5 has occurred 
in iteration, we replace to by y and go back to step 2. If 4>(to) > 5, we go to step 4. 

4. If 4>{to) > 5 and the inequality </>(to) < has not occured in iterations, we replace 
to by 3to and go back to step 2. If the inequality 4>{to) < has occured in some 
iteration before, we stop the iteration and use to as an initial guess in Newton's 
iterations for finding tg. 

4 Numerical experiments 

In this section results of some numerical experiments with ill-conditioned linear algebraic 
systems are reported. In all the experiments, by DSMG we denote the version of the DSM 
described in this paper, by VR we denote the Variational Regularization, implemented 
using the discrepancy principle, and by DSM- [2] we denote the method developed in [2j. 

4.1 A linear algebraic system for the computation of second derivatives 

Let us do some numerical experiments with linear algebraic systems arising in a numerical 
experiment of computing the second derivative of a noisy function. 

The problem is reduced to an integral equation of the first kind. A linear algebraic 
system is obtained by a discretization of the integral equation whose kernel K is Green's 
function 



< (j)(to) = \\e- toSs h s \\- 5 <5 



(44) 
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Here s,t S [0,1]. Using from [2J, we do some numerical experiments for solving 
from the linear algebraic system A^un = bj^g- I n the experiments the exact right-hand 
side is computed by the formula bjy = A^u^ when is given. In this test, un is 
computed by 

u N := (u(t Njl ),u(t Nj2 ), ....,u(t N>N )) , t Nji :=—, i = l,...,N, 

where u(t) is a given function. We use N = 10,20, 100 and 6/v,<5 = &/V + cat, where eAr 
is a random vector whose coordinates are independent, normally distributed, with mean 
and variance 1, and scaled so that [|ejv|| = <W||&iv||- This linear algebraic system is mildly 
ill-posed: the condition number of ^4ioo is 1-2158 x 10 4 . 

In Figure [H the difference between the exaction and solution obtained by the DSMG, 
VR and DSM-[2] are plotted. In these experiments, we used N = 100 and u(t) = sin(-7ri) 
with 5 re [ = 0.05 and 5 re [ = 0.01. Figure Q] shows that the results obtained by the VR and 
the DSM-[2] are very close to each other. The results obtained by the DSMG are much 
better than those by the DSM-[2] and by the VR. 




Figure 1: Plots of differences between the exact solution and solutions obtained by the 
DSMG, VR and DSM-ffl. 

Table [T] presents numerical results when N varies from 10 to 100, u(t) = sin(27ri), and 
t G [0, 1]. In this experiment the DSMG yields more accurate solutions than the DSM-[2] 
and the VR. The DSMG in this experiment takes more iterations than the DSM-[2] and 
the VR to get a solution. 

In this experiment the DSMG is implemented using the SVD of A obtained by the 
function svd in Matlab. As already mentioned, the SVD is a special case of the spectral 
decomposition (|3ip . It is expensive to compute the SVD, in general. However, there are 
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Table 1: Numerical results for computing second derivatives with 5 re i = 0.01. 







DSM 


DSM- 2] 




VR 


N 


^itcr 


\\us-y\ 2 


H-linsol 


ll««-»Ha 




\\ u $~y\ h 




llvlla 


II V II 2 


20 


9 


0.0973 


3 


0.1130 


6 


0.1079 


30 


5 


0.0831 


4 


0.1316 


6 


0.1160 


40 


7 


0.0488 


4 


0.1150 


6 


0.1045 


50 


9 


0.0614 


4 


0.1415 


6 


0.1063 


60 


6 


0.0419 


4 


0.0919 


6 


0.0817 


70 


9 


0.0513 


4 


0.0961 


6 


0.0842 


80 


6 


0.0418 


4 


0.1225 


6 


0.0981 


90 


7 


0.0287 


4 


0.0919 


7 


0.0840 


100 


7 


0.0248 


5 


0.0778 


7 


0.0553 



practically important problems where the spectral decomposition (13ip can be computed 
fast (see Section [4.21 below). These problems consist of deconvolution problems using the 
Fast Fourier Transform (FFTs). 

The conclusion from this experiment is: the DSMG may yield results with much better 
accuracy than the VR and DSM-[2]. Numerical experiments for various u(t) show that 
the DSMG competes favorably with the VR and the DSM- [2]. 

4.2 An application to image restoration 

The image degradation process can be modeled by the following equation: 

g s =g + w, g = h*f, \\w\\ < 5, (45) 

where h represents a convolution function that models the blurring that many imaging 
systems introduce. For example, camera defocus, motion blur, imperfections of the lenses, 
all these phenomenon can be modeled by choosing a suitable h. The functions g$, /, and 
w are the observed image, the original signal, and the noise, respectively. The noise w 
can be due to the electronics used (thermal and shot noise), the recording medium (film 
grain), or the imaging process (photon noise). 

In practice g, h and / in equation (|4"5l) are often given as functions of a discrete argu- 
ment and equation (|45p can be written in this case as 

oo 

gs,i = 9i + Wi= ^2 fj h i-j + w i' «GZ. (46) 

j=-oo 

Note that one (or both) signals fj and hj have compact support (finite length). Suppose 
that signal / is periodic with period N, i.e., fi+N = fi-, and hj = for j < and 
j > N. Assume that / is represented by a sequence /o, /j\r-l and h is represented 
by ho, fiN-i. Then the convolution h * f is periodic signal g with period N, and the 
elements of g are defined as 

N-l 

9i= ^2 h jf(i-j)modN, i = 0,l,...,N -I. (47) 

3=0 
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Here (i — j) modN is i—j modulo N. The discrete Fourier transform (DFT) of g is denned 
as the sequence 



N-l 

9k ■= W 

3=0 



-i2wjk/N 



k = 0,1, ...,N - 1. 



Denote g = (go, g~N-i) T - Then equation (|4"7|) implies 

9 = fh, fg := (Mo, Ml, fN-ihN-if ■ 

Let diag(a) denote a diagonal matrix whose diagonal is (do, ctjv-i] 
are zeros. Then equation (|48p can be rewritten as 

g = A/, A := diag(/t). 



(48) 

and other entries 
(49) 



Since ^4 is of the form (|3ip with U = V = I and 5 = diag(/i), one can use the DSMG 
method to solve equation (|4T)|) stably for /. 

The image restoration test problem we use is taken from [7]. This test problem was 
developed at the US Air Force Phillips Laboratory, Lasers and Imaging Directorate, Kirt- 
land Air Force Base, New Mexico. The original and blurred images have 256 x 256 pixels, 
and are shown in Figure [21 These data has been widely used in the literature for testing 
image restoration algorithms. 



Original 




Blured-Noisy image 



Figure 2: Original and Blurred-noisy images. 

Figure [3] plots the regularized images by the VR and the DSMG when 5 re i = 0.01. 
Again, with an input value for 5 re i, the observed blurred-noisy images is computed by 

98 = 9 + re i- n err, 

1 1 err 

where err is a vector with random entries normally distributed with mean and variance 
1. In this experiment, it took 5 and 8 iterations for the DSMG and the VR, respectively, 
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to yield numerical results. Prom Figure one concludes that the DSMG is comparable to 
the VR in terms of accuracy. The time of computation in this experiment is about the 
same for the VR and DSMG. 



VR DSM 





* 








Figure 3: Regularized images when noise level is 1%. 

Figure S] plots the regularized images by the VR and the DSMG when 5 re i = 0.05. It 
took 4 and 7 iterations for the DSMG and the VR, respectively, to yield numerical results. 
Figure [H shows that the images obtained by the DSMG and the VR are about the same. 

VR DSM 













Figure 4: Regularized images when noise level is 5%. 



The conclusions from this experiment are: the DSMG yields results with the same 
accuracy as the VR, and requires less iterations than the VR. The restored images by the 
DSM- [2] are about the same as those by the VR. 
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Remark 2 Equation (|45p can be reduced to equation ([48 p whenever one of the two func- 
tions / and h has compact support and the other is periodic. 

5 Concluding remarks 

A version of the Dynamical Systems Method for solving ill-conditioned linear algebraic 
systems is studied in this paper. An a priori and a posteriori stopping rules are formu- 
lated and justified. An algorithm for computing the solution in the case when a spectral 
decomposition of the matrix A is available is presented. Numerical results show that the 
DSMG, i.e., the DSM version developed in this paper, yields results comparable to those 
obtained by the VR and the DSM-[2] developed in [2], and the DSMG method may yield 
much more accurate results than the VR method. It is demonstrated in [7] that the rate 
of convergence of the Landweber method can be increased by using preconditioning tech- 
niques. The rate of convergence of the DSM version, presented in this paper, might be 
improved by a similar technique. The advantage of our method over the steepest descent 
in [7j is the following: the stopping time t$ can be found from a discrepancy principle by 
Newton's iterations for a wide range of initial guess to; when t$ is found one can compute 
the solution without any iterations. Also, our method requires less iterations than the 
steepest descent in [7J, which is an accelerated version of the Landweber method. 

References 

[1] Haircr, E., and N0rsett, S. P., and Wanner, G., Solving ordinary differential equation I, 
nonstiff problems, Springer, Berlin, 1987. 

[2] Hoang, N. S. and Ramm, A. G., Solving ill-conditioned linear algebraic systems by the dy- 
namical systems method (DSM), Inverse Problems in Sci. and Engineering, 16, N5, (2008), 
pp. 617 - 630. 

[3] Hoang, N. S. and Ramm, A. G, On stable numerical differentiation, Australian J. Math. 
Anal. Appl., 5, Nl, (2008), Article 5, pp.1-7. 

[4] Ivanov, V., Tanana, V., Vasin, V., Theory of ill-posed problems, VSP, Utrecht, 2002. 

[5] Lattes, J., Lions, J., Methode de quasi-reversibilite et applications, Dunod, Paris, 1967. 

[6] Morozov, Methods of solving incorrectly posed problems, Springer Verlag, New York, 1984. 

[7] Nagy, J. G. and Palmer, K. M., Steepest descent, CG, and iterative regularization of ill-posed 
problems, BIT Numerical Mathematics, 43(2003), 1003-1017. 

[8] Ramm, A. G. and Smirnova, A. B., On stable numerical differentiation, Mathem. of Compu- 
tation, 70, (2001), 1131-1153. 

[9] Ramm, A. G. and Smirnova, A. B., Stable numerical differentiation: when is it possible? 
Jour. Korean SIAM, 7, Nl, (2003), 47-61. 



16 



[10] Ramm, A. G. and Airapetyan, R., Dynamical systems and discrete methods for solving non- 
linear ill-posed problems, Appl. Math. Reviews, vol. 1, Ed. G. Anastassiou, World Sci. Pub- 
lishers, 2000, pp.491-536. 



[11] Ramm, A. G., Dynamical systems method for solving operator equations, Elsevier, Amster- 
dam, 2007. 

[12] Tautenhahn, U., On the asymptotical regularization of nonlinear ill-posed problems, Inverse 
Problems, 10 (1994) 1405-1418. 

[13] Vainikko, G., Veretennikov, A., Iterative processes in ill-posed problems, Nauka, Moscow, 
1996. 



17 



